home *** CD-ROM | disk | FTP | other *** search
/ Amiga Tools 3 / Amiga Tools 3.iso / grafik / raytracing / rayshade-4.0.6.3 / libray / libobj / blob.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-08-09  |  32.2 KB  |  1,074 lines
  1. /*
  2.  * blob.c
  3.  *
  4.  * Copyright (C) 1990, 1991, Mark Podlipec, Craig E. Kolb
  5.  * All rights reserved.
  6.  *
  7.  * This software may be freely copied, modified, and redistributed
  8.  * provided that this copyright notice is preserved on all copies.
  9.  *
  10.  * You may not distribute this software, in whole or in part, as part of
  11.  * any commercial product without the express consent of the authors.
  12.  *
  13.  * There is no warranty or other guarantee of fitness of this software
  14.  * for any purpose.  It is provided solely "as is".
  15.  *
  16.  * blob.c,v 4.1 1994/08/09 07:58:00 explorer Exp
  17.  *
  18.  * blob.c,v
  19.  * Revision 4.1  1994/08/09  07:58:00  explorer
  20.  * Bump version to 4.1
  21.  *
  22.  * Revision 1.1.1.1  1994/08/08  04:52:05  explorer
  23.  * Initial import.  This is a prerelease of 4.0.6enh3, or 4.1 possibly.
  24.  *
  25.  * Revision 2.1  1992/06/02  14:08:51  mcgregor
  26.  * Installed the individual metaball surface support.  A step up from the
  27.  * individual color specification support. !NB! Attenuation adjustment may
  28.  * look kludgy, but a clear model of how to "blend" attenuation values is
  29.  * not available right now.
  30.  *
  31.  * Revision 1.1.1.1  1992/05/08  17:50:59  mcgregor
  32.  * This baby is working!!!!!  I placed the metaball_surface_contrib
  33.  * call in BlobIntersect per Mark's suggestion.  Next-- SURFACES!
  34.  *
  35.  * Revision 1.2  1992/05/07  20:23:45  mcgregor
  36.  * This version works, again, at generating a multiple colored blob primitive
  37.  * with full texturing.  The bug pointed out by Mark Podlipec still exists:
  38.  * a shadowed (or is it non-shadowed?) surface uses the ambient value from
  39.  * the last trip through this code.  It has to do with BlobNormal _NOT_
  40.  * being called for some things.  I just want to save a working edition even
  41.  * if it is only 90% working.
  42.  *
  43.  * Revision 1.1  1992/04/27  20:17:23  mcgregor
  44.  * Initial revision
  45.  *
  46.  * Revision 4.0  91/07/17  14:36:07  kolb
  47.  * Initial version.
  48.  * 
  49.  */
  50. #include "libcommon/common.h"
  51. #include "libsurf/surface.h"
  52. #include "geom.h"
  53. #include "blob.h"
  54.  
  55. static Methods *iBlobMethods = NULL;
  56. static char blobName[] = "blob";
  57.  
  58. unsigned long BlobTests, BlobHits;
  59. Surface init_surf; /* Preserve the surface details as found in input file */
  60.  
  61. /*
  62.  * Blob/Metaball Description
  63.  *
  64.  * In this implementation a Blob is made up of a threshold T, and a 
  65.  * group of 1 or more Metaballs.  Each metaball 'i'  is defined by 
  66.  * its position (xi,yi,zi), its radius ri, and its strength si.
  67.  * Around each Metaball is a density function di(Ri) defined by:
  68.  *
  69.  *         di(Ri) =  c4i * Ri^4  +  c2i * Ri^2  + c0i     0 <= Ri <= ri
  70.  *         di(Ri) =  0                                    ri < Ri 
  71.  *
  72.  * where
  73.  *       c4i  = si / (ri^4)
  74.  *       c2i  = -(2 * si) / (ri^2)
  75.  *       c0i  = si
  76.  *       Ri^2 is the distance from a point (x,y,z) to the center of
  77.  *            the Metaball - (xi,yi,zi).
  78.  *
  79.  * The density function looks sort of like a W (Float-u) with the 
  80.  * center hump crossing the d-axis at  si  and the two humps on the
  81.  * side being tangent to the R-axis at  -ri  and  +ri. Only the 
  82.  * range  [0,ri] is being used. 
  83.  * I chose this function so that derivative = 0  at the points  0 and +ri. 
  84.  * This keeps the surface smooth when the densities are added. I also 
  85.  * wanted no  Ri^3  and  Ri  terms, since the equations are easier to 
  86.  * solve without them. This led to the symmetry about the d-axis and
  87.  * the derivative being equal to zero at -ri as well.
  88.  *
  89.  * The surface of the Blob is defined by:
  90.  *
  91.  *
  92.  *                  N
  93.  *                 ---  
  94.  *      F(x,y,z) = \    di(Ri)  = T
  95.  *                 /
  96.  *                 ---
  97.  *                 i=0
  98.  *
  99.  * where
  100.  *
  101.  *     di(Ri)    is   given above
  102.  *     Ri^2      =    (x-xi)^2  +  (y-yi)^2  + (z-zi)^2
  103.  *      N        is   number of Metaballs in Blob
  104.  *      T        is   the threshold
  105.  *    (xi,yi,zi) is   the center of Metaball i
  106.  *
  107.  */
  108.  
  109. /*****************************************************************************
  110.  * Create & return reference to a metaball.
  111.  * gnm 4-30-92  added reference to blob's surface
  112.  */
  113. Blob *
  114. BlobCreate(surf, T, mlist, npoints)
  115.      Surface *surf;               /* surface reference for the blob */
  116.      Float T;                     /* Treshold */
  117.      MetaList *mlist;             /* Pointer to linked list of metaballs */
  118.      int npoints;                 /* Number of metaballs in blob instance */
  119. {
  120.   Blob *blob;
  121.   int i;
  122.   MetaList *cur;
  123.   
  124.   /* 
  125.    * There has to be at least one metaball in the blob.
  126.    * Note: if there's only one metaball, the blob 
  127.    * will be a sphere.
  128.    */
  129.   if (npoints < 1)
  130.     {
  131.       RLerror(RL_WARN, "blob field not correct.\n");
  132.       return (Blob *)NULL;
  133.     };
  134.   
  135.   
  136.   blob = (Blob *)Malloc(sizeof(Blob));
  137.   
  138.   if((blob->surf = surf) != NULL)
  139.     {
  140.       /* save the original surface details ...er... oops this'll
  141.      only save the most *recent* blob instance's surface.
  142.      I have a plan to save each blob's (and metaballs' ??!)
  143.      surface, but I'll implement that later.
  144.      */
  145.       save_orig_master_surf(surf);
  146.     };
  147.   blob->T = T;
  148.   blob->list=(MetaVector *)Malloc((unsigned)(npoints*sizeof(MetaVector)));
  149.   blob->num = npoints;
  150.  
  151.   /* blob will, initially be single color, not multicolored. */
  152.   blob->MColor_Flag = FALSE;
  153.  
  154.   /*
  155.    * Initialize Metaball list
  156.    */
  157.   for(i=0;i<npoints;i++)
  158.     {
  159.       cur = mlist;
  160.       if (   ( (cur->mvec.c0 > -EPSILON) && (cur->mvec.c0 < EPSILON) )
  161.       || (cur->mvec.rs < EPSILON) )
  162.     {
  163.       RLerror(RL_WARN,"Degenerate metaball in blob (sr).\n");
  164.       return (Blob *)NULL;
  165.     };
  166.  
  167.       /* store radius squared */
  168.       blob->list[i].rs = cur->mvec.rs * cur->mvec.rs;
  169.  
  170.       /* Calculate and store coefficients for each metaball */
  171.       blob->list[i].c0 = cur->mvec.c0;
  172.       blob->list[i].c2 = -(2.0 * cur->mvec.c0) / blob->list[i].rs;
  173.       blob->list[i].c4 = cur->mvec.c0 /    (blob->list[i].rs * blob->list[i].rs);
  174.       blob->list[i].x = cur->mvec.x;
  175.       blob->list[i].y = cur->mvec.y;
  176.       blob->list[i].z = cur->mvec.z;
  177.  
  178.       blob->list[i].surf = cur->mvec.surf;
  179.  
  180.      /* store metaball black body color info
  181.       blob->list[i].mcolor.red = cur->mvec.mcolor.red;
  182.       blob->list[i].mcolor.green = cur->mvec.mcolor.green;
  183.       blob->list[i].mcolor.blue = cur->mvec.mcolor.blue;
  184.       */
  185.  
  186.       /* If any mcolor value is != 1.0 we need to assume that
  187.      this blob will be multi colored.
  188.      if((blob->list[i].mcolor.red != 1.0)
  189.      || (blob->list[i].mcolor.green != 1.0)
  190.      || (blob->list[i].mcolor.blue != 1.0))
  191.      */
  192.  
  193.       /* if the surf spec is non-NULL then we have a multisurfaced blob */
  194.       if(blob->list[i].surf != NULL)
  195.     blob->MColor_Flag = TRUE;
  196.       
  197.       mlist = mlist->next;
  198.       free((voidstar)cur);
  199.     };
  200.  
  201.   /* 
  202.    * Allocate room for Intersection Structures and
  203.    * Allocate room for an array of pointers to point to
  204.    * the Intersection Structures.
  205.    */
  206.   blob->ilist = (MetaInt *)Malloc( 2 * blob->num * sizeof(MetaInt));
  207.   blob->iarr = (MetaInt **)Malloc( 2 * blob->num * sizeof(MetaInt *));
  208.   return blob;
  209. }
  210.  
  211.  
  212.  
  213. /* :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
  214.  *
  215.  *  Function: save_orig_master_surf(surf)
  216.  *                Surface *surf;
  217.  *
  218.  *            Saves the original, master surface description(if any).
  219.  *
  220.  *  Errors:
  221.  *
  222.  *  Notes:  This needs some _SERIOUS_ reworking!
  223.  *
  224.  *  --------------------------------------------------------------------
  225.  *  who                  when                  what
  226.  *  --------------------------------------------------------------------
  227.  *  GN McGregor
  228.  *  --------------------------------------------------------------------
  229.  */
  230.  
  231. int
  232. save_orig_master_surf(surf)
  233.      Surface *surf;
  234. {
  235.   
  236.   
  237.   init_surf.name = surf->name;  /* Name */
  238.   init_surf.amb.r = surf->amb.r;  /* Ambient 'curve' */
  239.   init_surf.amb.g = surf->amb.g;  /* Ambient 'curve' */
  240.   init_surf.amb.b = surf->amb.b;  /* Ambient 'curve' */
  241.   init_surf.diff.r = surf->diff.r;  /* Diffuse reflection 'curve' */
  242.   init_surf.diff.g = surf->diff.g;  /* Diffuse reflection 'curve' */
  243.   init_surf.diff.b = surf->diff.b;  /* Diffuse reflection 'curve' */
  244.   init_surf.spec.r = surf->spec.r;  /* Specular reflection 'curve' */
  245.   init_surf.spec.g = surf->spec.g;  /* Specular reflection 'curve' */
  246.   init_surf.spec.b = surf->spec.b;  /* Specular reflection 'curve' */
  247.   init_surf.translu.r = surf->translu.r;  /* Diffuse transmission 'curve' */
  248.   init_surf.translu.g = surf->translu.g;  /* Diffuse transmission 'curve' */
  249.   init_surf.translu.b = surf->translu.b;  /* Diffuse transmission 'curve' */
  250.   init_surf.body.r = surf->body.r;  /* Specular transmission 'curve' */
  251.   init_surf.body.g = surf->body.g;  /* Specular transmission 'curve' */
  252.   init_surf.body.b = surf->body.b;  /* Specular transmission 'curve' */
  253.   init_surf.srexp = surf->srexp;  /* Specular reflection exponent */
  254.   init_surf.stexp = surf->stexp;  /* Specular transmission exponent */
  255.   init_surf.statten = surf->statten;  /* Specular transmission attenuation */
  256.   init_surf.index = surf->index;  /* Index of refraction */
  257.   init_surf.reflect = surf->reflect;  /* Specular reflectivity */
  258.   init_surf.transp = surf->transp;  /* Specular transmittance */
  259.   init_surf.translucency = surf->translucency;  /* Diffuse transmittance */ 
  260.   init_surf.noshadow = surf->noshadow;  /* No shadowing? */
  261.   init_surf.next = surf->next;  /* Next surface in list (if any) */
  262.  
  263.   return;
  264. }
  265.  
  266.  
  267. /* :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
  268.  *
  269.  *  Function: Methods *BlobMethods()
  270.  *
  271.  *
  272.  *  Errors:
  273.  *
  274.  *  Notes:
  275.  *
  276.  *  --------------------------------------------------------------------
  277.  *  who                  when                  what
  278.  *  --------------------------------------------------------------------
  279.  *  GN McGregor
  280.  *  --------------------------------------------------------------------
  281.  */
  282.  
  283. Methods *
  284. BlobMethods()
  285. {
  286.   if (iBlobMethods == (Methods *)NULL)
  287.     {
  288.       iBlobMethods = MethodsCreate();
  289.       iBlobMethods->create = (GeomCreateFunc *)BlobCreate;
  290.       iBlobMethods->methods = BlobMethods;
  291.       iBlobMethods->name = BlobName;
  292.       iBlobMethods->intersect = BlobIntersect;
  293.       iBlobMethods->normal = BlobNormal;
  294.       iBlobMethods->bounds = BlobBounds;
  295.       iBlobMethods->stats = BlobStats;
  296.       iBlobMethods->checkbounds = TRUE;
  297.       iBlobMethods->closed = TRUE;
  298.     };
  299.   return(iBlobMethods);
  300. }
  301.  
  302. /*****************************************************************************
  303.  * Function used by qsort() when sorting the Ray/Blob intersection list
  304.  */
  305. int
  306. MetaCompare(A,B)
  307.      char *A,*B;
  308. {
  309.   MetaInt **AA,**BB;
  310.   
  311.   AA = (MetaInt **) A;
  312.   BB = (MetaInt **) B;
  313.   if (AA[0]->bound == BB[0]->bound)
  314.     return(0);
  315.   if (AA[0]->bound <  BB[0]->bound)
  316.     return(-1);
  317.   return(1);  /* AA[0]->bound is > BB[0]->bound */
  318. }
  319.  
  320. /*****************************************************************************
  321.  * Ray/metaball intersection test.
  322.  */
  323.  
  324. int
  325. BlobIntersect(blob, ray, mindist, maxdist)
  326.      Blob *blob;
  327.      Ray *ray;
  328. Float mindist, *maxdist;
  329. {
  330.   double c[5], root;
  331.   Float dist;
  332.   MetaInt *ilist,**iarr;
  333.   register int i,j,inum;
  334.   extern void qsort();
  335.   int metaball_surface_contrib(); /* Modify surface details accord. to di(Ri) */
  336.   BlobTests++;
  337.   
  338.   ilist = blob->ilist;
  339.   iarr = blob->iarr;
  340.   
  341.   
  342.   /*
  343.    * The first step in calculating the Ray/Blob intersection is to 
  344.    * divide the Ray into intervals such that only a fixed set of 
  345.    * Metaballs contribute to the density function along that interval. 
  346.    * This is done by finding the set of intersections between the Ray 
  347.    * and each Metaball's Sphere/Region of influence, which has a 
  348.    * radius  ri  and is centered at (xi,yi,zi).
  349.    * Intersection information is kept track of in the MetaInt 
  350.    * structure and consists of:
  351.    *
  352.    *   type    indicates whether this intersection is the start(R_START)
  353.    *           of a Region or the end(R_END) of one. 
  354.    *   pnt     the Metaball of this intersection
  355.    *   bound   the distance from Ray origin to this intersection
  356.    *
  357.    * This list is then sorted by  bound  and used later to find the Ray/Blob 
  358.    * intersection.
  359.    */
  360.   
  361.   inum = 0;
  362.   for(i=0; i < blob->num; i++) /* Check all metaballs for intersection */
  363.     {
  364.       register Float xadj, yadj, zadj;
  365.       register Float b, t, rs;
  366.       register Float dmin,dmax;
  367.       
  368.       rs   = blob->list[i].rs;
  369.       xadj = blob->list[i].x - ray->pos.x;
  370.       yadj = blob->list[i].y - ray->pos.y;
  371.       zadj = blob->list[i].z - ray->pos.z;
  372.       
  373.       /*
  374.        * Ray/Sphere of Influence intersection
  375.        */
  376.       b = xadj * ray->dir.x + yadj * ray->dir.y + zadj * ray->dir.z;
  377.       t = b * b - xadj * xadj - yadj * yadj - zadj * zadj + rs;
  378.       
  379.       /* 
  380.        * Don't accept imaginary or single roots. A single root is a ray
  381.        * tangent to the Metaball's Sphere/Region. The Metaball's
  382.        * contribution to the overall density function at this point is
  383.        * zero anyway.
  384.        */
  385.       if (t > 0.0) /* we have two roots and an intersection*/
  386.     {
  387.       t = sqrt(t);
  388.       dmin = b - t;
  389.       /* 
  390.        * only interested in stuff in front of ray origin 
  391.        */
  392.             if (dmin < mindist + EPSILON) dmin = mindist + EPSILON;
  393.       dmax = b + t;
  394.       if (dmax > dmin) /* we have a valid Region */
  395.         {
  396.           /*
  397.            * Initialize min/start and max/end Intersections Structures
  398.            * for this Metaball
  399.            */
  400.           ilist[inum].type = R_START;
  401.           ilist[inum].pnt = i;
  402.           ilist[inum].bound = dmin;
  403.           for(j=0;j<5;j++) ilist[inum].c[j] = 0.0;
  404.           iarr[inum] = &(ilist[inum]);
  405.           inum++;
  406.           
  407.           ilist[inum].type = R_END;
  408.           ilist[inum].pnt = i;
  409.           ilist[inum].bound = dmax;
  410.           for(j=0;j<5;j++) ilist[inum].c[j] = 0.0;
  411.           iarr[inum] = &(ilist[inum]);
  412.           inum++;
  413.         }; /* End of valid Region */
  414.     }; /* End of two roots */
  415.     }; /* End of loop through metaballs */
  416.   
  417.   /*
  418.    * If there are no Ray/Metaball intersections there will 
  419.    * not be a Ray/Blob intersection. Exit now.
  420.    */
  421.   if (inum == 0)
  422.     {
  423.       return FALSE;
  424.     }
  425.   
  426.   /* 
  427.    * Sort Intersection list. No sense using qsort if there's only
  428.    * two intersections.
  429.    *
  430.    * Note: we actually aren't sorting the Intersection structures, but
  431.    * an array of pointers to the Intersection structures. 
  432.    * This is faster than sorting the Intersection structures themselves.
  433.    */
  434.   if (inum==2)
  435.     {
  436.       MetaInt *t;
  437.       if (ilist[0].bound > ilist[1].bound)
  438.     {
  439.       t = iarr[0];
  440.       iarr[0] = iarr[1];
  441.       iarr[1] = t;
  442.     };
  443.     }
  444.   else
  445.     qsort((voidstar)(iarr), (unsigned)inum, sizeof(MetaInt *), MetaCompare);
  446.   
  447.   /*
  448.    * Finding the Ray/Blob Intersection
  449.    * 
  450.    * The non-zero part of the density function for each Metaball is
  451.    * 
  452.    *   di(Ri) = c4i * Ri^4  +  c2i * Ri^2  +  c0i 
  453.    *
  454.    * To find find the Ray/Blob intersection for one metaball
  455.    * substitute for distance   
  456.    *
  457.    *     Ri^2 = (x-xi)^2 + (y-yi)^2 + (z-zi)^2
  458.    * 
  459.    * and then substitute the Ray equations:
  460.    *   
  461.    *     x  = x0 + x1 * t
  462.    *     y  = y0 + y1 * t
  463.    *     z  = z0 + z1 * t
  464.    *
  465.    * to get a big mess :^). Actually, it's a Quartic in t and it's fully 
  466.    * listed farther down. Here's a short version:
  467.    *
  468.    *   c[4] * t^4  +  c[3] * t^3  +  c[2] * t^2  +  c[1] * t  +  c[0]  =  T
  469.    *
  470.    * Don't forget that the Blob is defined by the density being equal to 
  471.    * the threshold T.
  472.    * To calculate the intersection of a Ray and two or more Metaballs, 
  473.    * the coefficients are calculated for each Metaball and then added 
  474.    * together. We can do this since we're working with polynomials.
  475.    * The points of intersection are the roots of the resultant equation.
  476.    *
  477.    * The algorithm loops through the intersection list, calculating
  478.    * the coefficients if an intersection is the start of a Region and 
  479.    * adding them to all intersections in that region.
  480.    * When it detects a valid interval, it calculates the 
  481.    * roots from the starting intersection's coefficients and if any of 
  482.    * the roots are in the interval, the smallest one is returned.
  483.    *
  484.    */
  485.   
  486.   { /* This curly brace is for variable scoping control. */
  487.     register Float *tmpc;
  488.     MetaInt *strt,*tmp;
  489.     register int istrt,itmp;
  490.     register int num,exitflag,inside;
  491.     int clear;
  492.     
  493.     /*
  494.      * Start the algorithm outside the first interval
  495.      */
  496.     inside = 0;
  497.     istrt = 0; 
  498.     strt = iarr[istrt];
  499.     if (strt->type != R_START)
  500.       RLerror(RL_WARN,"MetaInt sanity check FAILED!\n");
  501.     
  502.     /*
  503.      * Loop through intersection. If a root is found the code
  504.      * will return at that point.
  505.      */
  506.     while( istrt < inum ) /* inum ==> number of intersections aka iarr entries */
  507.       {
  508.     /* 
  509.      * Check for multiple intersections  at the same point.
  510.      * This is also where the coefficients are calculated
  511.      * and spread throughout that Metaball's sphere of
  512.      * influence.
  513.      */
  514.     do /* BEGIN: search for valid interval(s) */
  515.       {
  516.         /* find out which metaball */
  517.         i = strt->pnt;
  518.         /* only at starting boundaries do this */
  519.         if (strt->type == R_START)
  520.           {
  521.         register MetaVector *ml;
  522.         register Float a1,a0;
  523.         register Float xd,yd,zd;
  524.         register Float c4,c2,c0;
  525.         
  526.         /* multi colored blob variables */
  527.         register Float Ri2; /* Ri^2 */
  528.         register Float Ri4; /* Ri^4 */
  529.         register Float DF;  /* di(Ri) / T */
  530.         
  531.         /* we're inside */
  532.         inside++;
  533.         
  534.         /*  As promised, the full equations
  535.          *
  536.          *   c[4] = c4*a2*a2; 
  537.          *   c[3] = 4.0*c4*a1*a2;
  538.          *   c[2] = 4.0*c4*a1*a1 + 2.0*c4*a2*a0 + c2*a2;
  539.          *   c[1] = 4.0*c4*a1*a0 + 2.0*c2*a1;
  540.          *   c[0] = c4*a0*a0 + c2*a0 + c0;
  541.          *
  542.          * where
  543.          *        a2 = (x1*x1 + y1*y1 + z1*z1) = 1.0 because the ray 
  544.          *                                           is normalized
  545.          *        a1 = (xd*x1 + yd*y1 + zd*z1)
  546.          *        a0 = (xd*xd + yd*yd + zd*zd)
  547.          *        xd = (x0 - xi)
  548.          *        yd = (y0 - yi)
  549.          *        zd = (z0 - zi)
  550.          *        (xi,yi,zi) is center of Metaball
  551.          *        (x0,y0,z0) is Ray origin
  552.          *        (x1,y1,z1) is normalized Ray direction
  553.          *        c4,c2,c0   are the coefficients for the
  554.          *                       Metaball's density function
  555.          *
  556.          */
  557.         /*
  558.          * Some compilers would recalculate
  559.          * this each time its used.
  560.          */
  561.         ml = &(blob->list[i]);
  562.         
  563.         xd = ray->pos.x - ml->x;
  564.         yd = ray->pos.y - ml->y;
  565.         zd = ray->pos.z - ml->z;
  566.         a1 = (xd * ray->dir.x + yd * ray->dir.y + zd * ray->dir.z);
  567.         a0 = (xd * xd + yd * yd + zd * zd);
  568.         c0 = ml->c0;
  569.         c2 = ml->c2;
  570.         c4 = ml->c4;
  571.         c[4] = c4;
  572.         c[3] = 4.0*c4*a1;
  573.         c[2] = 2.0*c4*(2.0*a1*a1 + a0) + c2;
  574.         c[1] = 2.0*a1*(2.0*c4*a0 + c2);
  575.         c[0] = c4*a0*a0 + c2*a0 + c0;
  576.                 
  577.         /* 
  578.          * Since we've gone through the trouble of calculating the 
  579.          * coefficients, put them where they'll be used.
  580.          * Starting at the current intersection(It's also the start of
  581.          * a region), add the coefficients to each intersection until
  582.          * we reach the end of this region.
  583.          */
  584.         itmp = istrt; 
  585.         tmp = strt;
  586.         while( (tmp->pnt != i) || (tmp->type != R_END) )
  587.           {
  588.             tmpc = tmp->c;
  589.             for(j=0;j<5;j++)
  590.               *tmpc++ += c[j];
  591.             itmp++; 
  592.             tmp = iarr[itmp];
  593.           }
  594.         
  595.           } /* End of start of a Region */ 
  596.         /* 
  597.          * Since the intersection wasn't the start
  598.          * of a region, it must the end of one.
  599.          */
  600.         else
  601.           inside--;
  602.         
  603.         /*
  604.          * If we are inside a region(or regions) and the next
  605.          * intersection is not at the same place, then we have
  606.          * the start of an interval. Set the exitflag.
  607.          */
  608.         if ((inside > 0)
  609.         && (strt->bound != iarr[istrt+1]->bound) )
  610.           exitflag = 1;
  611.         else 
  612.           /* move to next intersection */
  613.           {
  614.         istrt++; 
  615.         strt = iarr[istrt];
  616.         exitflag = 0;
  617.           }
  618.         /* 
  619.          * Check to see if we're at the end. If so then exit with
  620.          * no intersection found.
  621.          */
  622.         if (istrt >= inum)
  623.           {
  624.         return FALSE;
  625.           }
  626.       }
  627.     while(!exitflag); /* END: of Search for valid interval (do-while) */
  628.     
  629.     /*
  630.      * Find Roots along this interval
  631.      */
  632.     
  633.     /* get coefficients from Intersection structure */
  634.     tmpc = strt->c;
  635.     for(j=0;j<5;j++)
  636.       c[j] = *tmpc++;
  637.     
  638.     /* Don't forget to put in threshold */
  639.     c[0] -= blob->T;
  640.     
  641.     /* Use Graphic Gem's Quartic Root Finder */
  642.     num = FindRoot(c,4,strt->bound,iarr[istrt+1]->bound,&root,-1);
  643.     
  644.  
  645.             if ( num && (root > mindist) && (root < *maxdist) ) 
  646.             {
  647.                  *maxdist = root;
  648.                  BlobHits++;
  649.                 if(blob->MColor_Flag)
  650.                     metaball_surface_contrib(blob,ray,root);
  651.                  return TRUE;
  652.             }
  653.     
  654.     /* 
  655.      * Move to next intersection
  656.      */
  657.     istrt++;
  658.     strt = iarr[istrt];
  659.     
  660.       }; /* End of loop through Intersection List */
  661.   }; /* End of Solving for Ray/Blob Intersection */
  662.   /* 
  663.    * return negative
  664.    */
  665.   return FALSE;
  666. }
  667.  
  668.  
  669. /***********************************************
  670.  * Find the Normal of a Blob at a given point
  671.  *
  672.  * The normal of a surface at a point (x0,y0,z0) is the gradient of that
  673.  * surface at that point. The gradient is the vector 
  674.  *
  675.  *            Fx(x0,y0,z0) , Fy(x0,y0,z0) , Fz(x0,y0,z0)
  676.  *
  677.  * where Fx(),Fy(),Fz() are the partial dervitives of F() with respect
  678.  * to x,y and z respectively. Since the surface of a Blob is made up
  679.  * of the sum of one or more polynomials, the normal of a Blob at a point
  680.  * is the sum of the gradients of the individual polynomials at that point.
  681.  * The individual polynomials in this case are di(Ri) i = 0,...,num .
  682.  *
  683.  * To find the gradient of the contributing polynomials
  684.  * take di(Ri) and substitute U = Ri^2;
  685.  *
  686.  *      di(U)    = c4i * U^2  +  c2i * U  + c0i
  687.  *
  688.  *      dix(U)   = 2 * c4i * Ux * U  +  c2i * Ux
  689.  *
  690.  *        U      = (x-xi)^2 + (y-yi)^2 + (z-zi)^2 
  691.  *
  692.  *        Ux     = 2 * (x-xi) 
  693.  *
  694.  * Rearranging we get
  695.  *  
  696.  *    dix(x0,y0,z0) = 4 * c4i * xdi * disti + 2 * c2i * xdi
  697.  *    diy(x0,y0,z0) = 4 * c4i * ydi * disti + 2 * c2i * ydi
  698.  *    diz(x0,y0,z0) = 4 * c4i * zdi * disti + 2 * c2i * zdi
  699.  * 
  700.  * where
  701.  *         xdi       =   x0-xi
  702.  *         ydi       =   y0-yi
  703.  *         zdi       =   y0-yi
  704.  *        disti      =   xdi*xdi + ydi*ydi + zdi*zdi   
  705.  *       (xi,yi,zi)  is  the center of Metaball i
  706.  *       c4i,c2i,c0i are the coefficients of Metaball i's density function
  707.  */
  708.  
  709. int
  710. BlobNormal(blob, pos, nrm, gnrm)
  711.      Blob *blob;
  712.      Vector *pos, *nrm, *gnrm;
  713. {
  714.   register int i;
  715.   int clear;
  716.   
  717.   /*  
  718.    * Initialize normals to zero 
  719.    */
  720.   nrm->x = nrm->y = nrm->z = 0.0;
  721.   
  722.   /*
  723.    * Loop through Metaballs. If the point is within a Metaball's
  724.    * Sphere of influence, calculate the gradient and add it to the
  725.    * normals
  726.    */
  727.   for(i=0;i < blob->num; i++)
  728.     {
  729.       register MetaVector *sl;
  730.       register Float dist,xd,yd,zd;
  731.       
  732.       sl = &(blob->list[i]);
  733.       xd = pos->x - sl->x;
  734.       yd = pos->y - sl->y;
  735.       zd = pos->z - sl->z;
  736.       
  737.       dist  = xd*xd + yd*yd + zd*zd;
  738.       if (dist <= sl->rs )
  739.     {
  740.       register Float temp;
  741.       register Float Ri4, Ri2, DF;
  742.       
  743.       /* temp is negative so normal points out of blob */
  744.       temp = -2.0 * (2.0 * sl->c4 * dist  +  sl->c2);
  745.       nrm->x += xd * temp;
  746.       nrm->y += yd * temp;
  747.       nrm->z += zd * temp;
  748.       
  749.     };
  750.     };
  751.   (void)VecNormalize(nrm);
  752.   *gnrm = *nrm;
  753.   return FALSE;
  754. }
  755.  
  756.  
  757. /*****************************************************************************
  758.  * Calculate the extent of the Blob
  759.  */
  760. void
  761. BlobBounds(blob, bounds)
  762.      Blob *blob;
  763.      Float bounds[2][3];
  764. {
  765.   Float r;
  766.   Float minx,miny,minz,maxx,maxy,maxz;
  767.   int i;
  768.   
  769.   /*
  770.    * Loop through Metaballs to find the minimun and maximum in each
  771.    * direction.
  772.    */
  773.   for( i=0;  i< blob->num;  i++)
  774.     {
  775.       register Float d;
  776.       
  777.       r = sqrt(blob->list[i].rs);
  778.       if (i==0)
  779.     {
  780.       minx = blob->list[i].x - r;
  781.       miny = blob->list[i].y - r;
  782.       minz = blob->list[i].z - r;
  783.       maxx = blob->list[i].x + r;
  784.       maxy = blob->list[i].y + r;
  785.       maxz = blob->list[i].z + r;
  786.     }
  787.       else
  788.     {
  789.       d = blob->list[i].x - r; 
  790.       if (d<minx) minx = d;
  791.       d = blob->list[i].y - r; 
  792.       if (d<miny) miny = d;
  793.       d = blob->list[i].z - r; 
  794.       if (d<minz) minz = d;
  795.       d = blob->list[i].x + r; 
  796.       if (d>maxx) maxx = d;
  797.       d = blob->list[i].y + r; 
  798.       if (d>maxy) maxy = d;
  799.       d = blob->list[i].z + r; 
  800.       if (d>maxz) maxz = d;
  801.     };
  802.     };
  803.   bounds[LOW][X]  = minx;
  804.   bounds[HIGH][X] = maxx;
  805.   bounds[LOW][Y]  = miny;
  806.   bounds[HIGH][Y] = maxy;
  807.   bounds[LOW][Z]  = minz;
  808.   bounds[HIGH][Z] = maxz;
  809.   return;
  810. }
  811.  
  812. char *
  813. BlobName()
  814. {
  815.   return blobName;
  816. }
  817.  
  818. void
  819. BlobStats(tests, hits)
  820.      unsigned long *tests, *hits;
  821. {
  822.   *tests = BlobTests;
  823.   *hits = BlobHits;
  824.   return;
  825. }
  826.  
  827. void
  828. BlobMethodRegister(meth)
  829.      UserMethodType meth;
  830. {
  831.   if (iBlobMethods)
  832.     iBlobMethods->user = meth;
  833.   return;
  834. }
  835.  
  836.  
  837.  
  838.  
  839. /* :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
  840.  *
  841.  *  Function: restore_orig_surf_descrip(blob) 
  842.  *            restores the original surface details as they were passed
  843.  *            into Rayshade for the blob instance.
  844.  *  Errors:  None at this time.
  845.  *
  846.  *  Notes: This _SHOULD_ bullet proof the code, but where?????
  847.  *
  848.  *  --------------------------------------------------------------------
  849.  *  who                  when                  what
  850.  *  --------------------------------------------------------------------
  851.  *  GN McGregor       6-May-92           development revisions
  852.  *  --------------------------------------------------------------------
  853.  */
  854.  
  855. int
  856. restore_orig_surf_descrip(blob)
  857.      Blob *blob;
  858. {
  859.   Surface *surf;
  860.   /* As the routine name suggests we are restoring the original surface
  861.      description as it was handed in to the program for the blob.  While
  862.      this may be a bit of overkill, it does make this part of the code a
  863.      bit more bullet proof.
  864.      */
  865.   
  866.   if((surf = blob->surf))
  867.     {
  868.       surf->name = init_surf.name;  /* Name */
  869.       surf->amb.r = init_surf.amb.r;  /* Ambient 'curve' */
  870.       surf->amb.g = init_surf.amb.g;  /* Ambient 'curve' */
  871.       surf->amb.b = init_surf.amb.b;  /* Ambient 'curve' */
  872.       surf->diff.r = init_surf.diff.r;  /* Diffuse reflection 'curve' */
  873.       surf->diff.g = init_surf.diff.g;  /* Diffuse reflection 'curve' */
  874.       surf->diff.b = init_surf.diff.b;  /* Diffuse reflection 'curve' */
  875.       surf->spec.r = init_surf.spec.r;  /* Specular reflection 'curve' */
  876.       surf->spec.g = init_surf.spec.g;  /* Specular reflection 'curve' */
  877.       surf->spec.b = init_surf.spec.b;  /* Specular reflection 'curve' */
  878.       surf->translu.r = init_surf.translu.r;  /* Diffuse transmission 'curve' */
  879.       surf->translu.g = init_surf.translu.g;  /* Diffuse transmission 'curve' */
  880.       surf->translu.b = init_surf.translu.b;  /* Diffuse transmission 'curve' */
  881.       surf->body.r = init_surf.body.r;  /* Specular transmission 'curve' */
  882.       surf->body.g = init_surf.body.g;  /* Specular transmission 'curve' */
  883.       surf->body.b = init_surf.body.b;  /* Specular transmission 'curve' */
  884.       surf->srexp = init_surf.srexp;  /* Specular reflection exponent */
  885.       surf->stexp = init_surf.stexp;  /* Specular transmission exponent */
  886.       surf->statten = init_surf.statten;  /* Specular transmission attenuation */
  887.       surf->index = init_surf.index;  /* Index of refraction */
  888.       surf->reflect = init_surf.reflect;  /* Specular reflectivity */
  889.       surf->transp = init_surf.transp;  /* Specular transmittance */
  890.       surf->translucency = init_surf.translucency;  /* Diffuse transmittance */ 
  891.       surf->noshadow = init_surf.noshadow;  /* No shadowing? */
  892.       surf->next = init_surf.next;  /* Next surface in list (if any) */
  893.       
  894.     };
  895.   return(GOOD);
  896. }
  897.  
  898.  
  899. /* :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
  900.  *
  901.  *  Function: metaball_surface_contrib(blob, ray, dist)
  902.  *
  903.  *
  904.  *  Errors:
  905.  *
  906.  *  Notes:
  907.  *
  908.  *  --------------------------------------------------------------------
  909.  *  who                  when                  what
  910.  *  --------------------------------------------------------------------
  911.  *  GN McGregor
  912.  *  --------------------------------------------------------------------
  913.  */
  914. int
  915. metaball_surface_contrib(blob, ray, dist)
  916.      Blob *blob;
  917.      Ray *ray;
  918.      Float dist;
  919. {
  920.   Vector intersect_pos;
  921.   int i;
  922.   
  923.   /* We have the distance from the ray's origin to the point of intersection.
  924.      We have to find the x, y, z values for the point equal to dist
  925.      along the ray.The intersection x, y, z values can be drived by 
  926.      extending the unit vector like so...
  927.      */
  928.  
  929.   VecAddScaled(ray->pos, dist, ray->dir, &intersect_pos);
  930.  
  931.   /* Prepare to modify the surface */
  932.   tweak_blob_surf(1,  blob->surf, 1.0, &(blob->list[0]));
  933.  
  934.   for(i=0;i < blob->num; i++) /* Scan metaballs for surface contributions */
  935.     {
  936.       register MetaVector *sl;
  937.       register Float range,xd,yd,zd;
  938.       register Float Ri2;  /* R(i)^2 in Mark's Equation up above. */
  939.       
  940.       sl = &(blob->list[i]);
  941.       xd = intersect_pos.x - sl->x;
  942.       yd = intersect_pos.y - sl->y;
  943.       zd = intersect_pos.z - sl->z;
  944.       
  945.       Ri2  = xd*xd + yd*yd + zd*zd;
  946.       if (Ri2 <= sl->rs )
  947.     {
  948.       register Float DF;
  949.  
  950.       DF = (((sl->c4*Ri2) + sl->c2)*Ri2 + sl->c0) / blob->T;
  951.       /*DF = ((sl->c4*Ri4) + (sl->c2*Ri2) + sl->c0) / blob->T;*/
  952.       tweak_blob_surf(0, blob->surf, DF, sl->surf);
  953.     }; /* END metaball density calc. */
  954.     };
  955.   
  956.   return(GOOD);
  957. }
  958.  
  959.  
  960.  
  961. /* :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
  962.  *
  963.  *  Function: tweak_blob_surf(clear, surf, df, mcolor)
  964.  *
  965.  *  This routine will fiddle with those components of the surface that
  966.  *    have RGB values as part of their specs.  The more sophisticated
  967.  *    concept and operation would involve allowing each metaball to possess
  968.  *    its own surface.  After this exercise that possibility may prove easier
  969.  *    and more useful.
  970.  *
  971.  *
  972.  *  Errors:
  973.  *
  974.  *  Notes:
  975.  *
  976.  *  --------------------------------------------------------------------
  977.  *  who                  when                  what
  978.  *  --------------------------------------------------------------------
  979.  *  GN McGregor     6-May-92          development revisions
  980.  *  --------------------------------------------------------------------
  981.  */
  982.  
  983. int
  984. tweak_blob_surf(clear, surf, df, msurf)
  985.      int clear;       /* flag-signal to init surface spec prior to tweaking */
  986.      Surface *surf;   /* pointer to surface spec. */
  987.      Float df;        /* density fraction-> contribution weight of metaball */
  988.      Surface *msurf;  /* metaball's surface spec. */
  989. {
  990.   
  991.   if(surf)
  992.     {
  993.       if(clear) /* We want to prepare the surface for its dynamic update. */
  994.     {
  995.       surf->amb.r = 0.0;
  996.       surf->amb.g = 0.0;
  997.       surf->amb.b = 0.0;
  998.       
  999.       surf->diff.r = 0.0;
  1000.       surf->diff.g = 0.0;
  1001.       surf->diff.b = 0.0;
  1002.       
  1003.       surf->spec.r = 0.0;
  1004.       surf->spec.g = 0.0;
  1005.       surf->spec.b = 0.0;
  1006.       
  1007.       surf->translu.r = 0.0;
  1008.       surf->translu.g = 0.0;
  1009.       surf->translu.b = 0.0;
  1010.       
  1011.       surf->body.r = 0.0;
  1012.       surf->body.g = 0.0;
  1013.       surf->body.b = 0.0;
  1014.  
  1015.       surf->srexp = 0.0;
  1016.       surf->stexp = 0.0;
  1017.       surf->statten = 0.0;
  1018.       surf->index = 0.0;
  1019.       surf->reflect = 0.0;
  1020.       surf->transp = 0.0;
  1021.       surf->translucency = 0.0;
  1022.     }
  1023.       else    /* We are updating the surface components. */
  1024.     {
  1025.       /* The new surface component color is the weighted sum of 
  1026.          each of the contributing  metaballs colors.  The  weighting
  1027.          factor is fractiobnal contribution of the metaball to the
  1028.          blob's density function. ( i.e. d(i) / T ).
  1029.          */
  1030.       /* Trying out: use surface specs as a scale value. */
  1031.       
  1032.       surf->amb.r += (msurf->amb.r * df) * init_surf.amb.r;
  1033.       surf->amb.g += (msurf->amb.g * df) * init_surf.amb.g;
  1034.       surf->amb.b += (msurf->amb.b * df) * init_surf.amb.b;
  1035.       
  1036.       surf->diff.r += (msurf->diff.r * df) * init_surf.diff.r;
  1037.       surf->diff.g += (msurf->diff.g * df) * init_surf.diff.g;
  1038.       surf->diff.b += (msurf->diff.b * df) * init_surf.diff.b;
  1039.       
  1040.       surf->spec.r += (msurf->spec.r * df) * init_surf.spec.r;
  1041.       surf->spec.g += (msurf->spec.g * df) * init_surf.spec.g;
  1042.       surf->spec.b += (msurf->spec.b * df) * init_surf.spec.b;
  1043.       
  1044.       surf->translu.r += (msurf->translu.r * df) * init_surf.translu.r;
  1045.       surf->translu.g += (msurf->translu.g * df) * init_surf.translu.g;
  1046.       surf->translu.b += (msurf->translu.b * df) * init_surf.translu.b;
  1047.       
  1048.       surf->body.r += (msurf->body.r * df) * init_surf.body.r;
  1049.       surf->body.g += (msurf->body.g * df) * init_surf.body.g;
  1050.       surf->body.b += (msurf->body.b * df) * init_surf.body.b;
  1051.  
  1052.       /* Model physical behavior as sum of weighted contributions from
  1053.          metaballs.
  1054.          */
  1055.       surf->srexp += (msurf->srexp * df);/* + init_surf.srexp;*/
  1056.       surf->stexp += (msurf->stexp * df);/* + init_surf.stexp;*/
  1057.  
  1058.       /* Attenuation _is_ a problem... Let's try a method where
  1059.          the attenuation is equal to the maximum value we see.
  1060.        */
  1061.       if(msurf->statten > init_surf.statten)
  1062.         surf->statten = (msurf->statten > surf->statten) ? msurf->statten : surf->statten;
  1063.       else
  1064.         surf->statten = (init_surf.statten > surf->statten) ? init_surf.statten : surf->statten;
  1065.  
  1066.       surf->index += (msurf->index * df);/* * init_surf.index;*/
  1067.       surf->reflect += (msurf->reflect * df);/* * init_surf.reflect;*/
  1068.       surf->transp += (msurf->transp * df);/* * init_surf.transp;*/
  1069.       surf->translucency += (msurf->translucency * df);/* * init_surf.translucency;*/
  1070.     };
  1071.     };
  1072.   return(GOOD);
  1073. }
  1074.